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We present an algorithm to construct meshes suitable for space-time discontinuous Galerkin finite-element methods. 
Our method generalizes and improves the 'Tent Pitcher' algorithm of Ungor and Sheffer. Given an arbitrary 
simplicially meshed domain X of any dimension and a time interval [0,T], our algorithm builds a simplicial mesh of 
the space-time domain X x [0,T], in constant time per element. Our algorithm avoids the limitations of previous 
methods by carefully adapting the durations of space-time elements to the local quality and feature size of the 
underlying space mesh. 
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1. INTRODUCTION 

Many simulation problems consider the behavior of an 
object or region of space over time. The most common 
finite element methods for this class of problem use 
a meshing procedure to discretize space, yielding a 
system of ordinary differential equations in time. A 
time-marching or time-integration scheme is then used 
to advance the solution over a series of fixed time steps. 
In general, a distinct spatial mesh may be required at 
each time step, due to the requirements of an adaptive 
analysis scheme or to track a moving boundary or 
interface within the domain. 

A relatively new approach to such simulations sug- 
gests directly meshing in space-time |)J |l^, F^2[. For 
example, a four-dimensional space-time mesh would 
be required to simulate an evolving three-dimensional 
domain. Usually, the time dimension is not treated 
in the same way as the spatial dimensions, in part 



because it can be scaled independently. Moreover, 
the numerical methods that motivate our research 
impose additional geometric constraints on the meshes 
to support a linear-time solution strategy. Thus, 
traditional meshing techniques do not apply. 

In this paper, we develop the first algorithm to build 
graded space-time meshes over arbitrary simplicially 
meshed domains in arbitrary dimensions. Our algo- 
rithm does not impose a fixed global time step on 
the mesh; rather, the duration of each space-time 
element depends on the local feature size and quality 
of the underlying space mesh. Our approach is a 
generalization of the 'Tent Pitcher' algorithm of Ungor 
and Sheffer j^], but avoids the restrictions of that 
method by imposing some additional constraints. Our 
algorithm builds space-time meshes in constant time 
per element. 

The paper is organized as follows. In Section FJ, 
we formalize the space-time meshing problem and 
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Figure 1. A space-time discontinuous Galerkin finite element mesh. 



describe several previous results. Section g explains 
the high-level advancing front strategy of our meshing 
algorithm. In Sections ^ and ^| we develop our 
algorithm for building three-dimensional space-time 
meshes over triangulated planar domains. We general- 
ize our algorithm to higher dimensions in Section [| In 
Section M, we describe our implementation and present 
some experimental results. Finally, we conclude in 
Section H by suggesting several directions for further 
research. 

2. SPACE-TIME DISCONTINUOUS 
GALERKIN MESHING 

The formulation of our space-time meshing problem 
relies on the notions of domain of influence and domain 
of dependence. Imagine dropping a pebble into a 
pond; over time, circular waves expand outward from 
the point of impact. These waves sweep out a cone in 
space-time, called the domain of influence of the event. 

More generally, we say that a point p in space-time de- 
pends on another point q if the salient physical param- 
eters at p (temperature, pressure, stress, momentum, 
etc.) can depend on the corresponding parameters at q, 
that is, if changing the conditions at q could change the 
conditions at p. The domain of influence of p is the set 
of points that depend on p; symmetrically, the domain 
of dependence is the set of points that p depends 
on. At least infinitesimally, these domains can be 
approximated by a pair of circular cones with common 
apex p. For isotropic problems without material flow, 
this double cone can described by a scalar wave speed 
c(p) € IR, which specifies how quickly the radius of the 



cones grows as a function of time. If the characteristic 
equations of the analysis problem are linear and the 
material properties are homogeneous, the wave speed 
is constant throughout the entire space-time domain; 
in this case, we can choose an appropriate time scale so 
that c(p) — 1 everywhere. For more general problems, 
the wave speed varies across space-time as a function 
of other physical parameters, and may even be part of 
the numerical solution. 

These notions extend to finite element meshes in 
space-time. We say that an element A in space-time 
depends on another element A' if any point p 6 A 
depends on any point q G A'. This relation naturally 
defines a directed dependency graph whose vertices are 
the elements of the mesh. Two elements in the mesh 
are coupled if they lie on a common directed cycle in 
(the transitive closure of) the dependency graph. 

Space-time discontinuous Galerkin (DG) methods 
have been proposed by Richter |12| , Lowrie et al. 
and Yin et al. J2^] for solving systems of nonlinear 
hyperbolic partial differential equations. These 
methods provide a linear-time element-by-element 
solution, avoiding the need to solve a large system of 
equations, provided no two elements in the underlying 
space-time mesh are coupled. In particular, every 
pair of adjacent elements must satisfy the so-called 
cone constraint: Any boundary facet between two 
neighboring elements separates the cone of influence 
from the cone of dependence of any point on the 
facet. See Figure |^. Intuitively, if a boundary facet 
satisfies the cone constraint, information can only 
flow in one direction across that facet. In a totally 
decoupled mesh, the dependency graph describes a 
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Figure 2. The cone constraint: Any boundary facet separates 
the domain of influence (above) from the domain of depen- 
dence (below). 

partial order on the elements, and the solution can 
be computed by considering the elements one at a 
time according to any linear extension of this partial 
order. Alternatively, the solutions within any set of 
incomparable elements can be computed in parallel. 

Discontinuous Galerkin methods impose no a priori 
restrictions on the shape of the individual elements; 
mixed meshes with tetrahedral, hexahedral, pyrami- 
dal, and other element shapes are acceptable. How- 
ever, it is usually more convenient to work with very 
simple convex elements such as simplices. Experience 
indicates that ill-conditioning is likely if the elements 
are non-convex, and subdividing non-convex regions 
into simple convex elements is useful for efficient 
integration. (For further background on DG methods, 
we refer the reader to the recent book edited by 
Cockburn, Karniadakis, and Shu |J, which contains 
both a general survey H and several papers describing 
space-time DG methods and their applications.) 

To construct an efficient mesh with convex elements, 
we have found it preferable to relax the cone con- 
straint in the following way. We construct a mesh 
of simplicial elements, but not all facets meet the cone 
constraint. Instead, elements are grouped into patches 
(of bounded size). The boundary facets between 
patches by definition satisfy the cone constraint, so 
patches are partially ordered by dependence, and can 
be solved independently. 

However, the internal facets between simplicial ele- 
ments within a patch may violate the cone constraint. 
Thus, DG methods require the elements within the 
patch to be solved simultaneously. Since each patch 
contains a constant number of elements, the system 
of equations within it has constant size, which implies 
that we can still solve the underlying numerical prob- 
lem in linear time by considering the patches one at a 
time. 

Richter jl^] observed that the dissipation of DG 
methods increases as the slope of boundary facets 
decreases below the local wave speed. Thus, our goal is 



to construct an efficient simplicial mesh, grouped into 
patches each containing few simplices, such that the 
boundary facets of each patch are as close as possible 
to the cone constraint without violating it. 

Previous Results 

Most previous space-time meshing algorithms con- 
struct a single mesh layer between two space-parallel 
planes and repeat this layer (or its reflection) at 
regular intervals to fill the simulation domain. The 
exact construction method depends on the type of un- 
derlying space mesh. For example, given a structured 
quad space mesh, the space-time meshing algorithm 
of Lowrie et al. M constructs a layer of pyramids and 
tetrahedra. Similarly, Ungor et al. j^, build a 
single layer of tetrahedra and pyramids over an acute 
triangular mesh, and Sheffer et al. jl4|, |2l] describe 
an algorithm to build a single layer of hexahedra over 
any (unstructured) quad mesh. All such layer-based 
approaches suffer from a global time step imposed by 
the smallest element in the underlying space mesh. 
This requirement increases the number of elements 
in the mesh, making the DG method less efficient; 
it also increases the numerical error of the solution, 
since many internal facets must lie significantly below 
the constraint cone. 

A few recent algorithms do not impose a global time 
step, but instead allows the durations of space-time 
elements to depend on the size of the underlying 
elements of the ground mesh. The first such algo- 
rithm, due to Ungor et al. pc| ], builds a triangular 
mesh for a (1 + 1) -dimensional space-time domain by 
intersecting the constraint cones at neighboring nodes. 
This method does not easily generalize to higher 
dimensions. The most general space-time meshing 
algorithm to date is the 'Tent Pitcher' algorithm of 
Ungor and Sheffer . Given a simplicial space mesh 
in any fixed dimension, where every dihedral angle is 
strictly less than 90° , Tent Pitcher constructs a space- 
time mesh of arbitrary duration. Moreover, if every 
dihedral angle in the space mesh is larger than some 
positive constant, each patch in the space-time mesh 
consists of a constant number of simplices. 

Unfortunately, the acute simplicial meshes that Tent 
Pitcher requires are difficult to construct, if not 
impossible, except in a few special cases. Bern 
et al. H describe two methods for building an acute 
triangular mesh for an arbitrary planar point set, and 
methods are known for special planar domains such 
as triangles [jilt , squares Q and some classes of 
polygons l^, |10| . However, no method is known for 
general planar domains or even for point sets in higher 
dimensions. It is an open problem whether the cube 
has an acute triangulation; see Jl7| for recent related 
results. 
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New Results 

In this paper, we present a generalization of the Tent 
Pitcher algorithm that extends any simplicial space 
mesh in TR d , for any d > 1, into a space-time mesh of 
arbitrary duration. Like the Tent Pitcher algorithm, 
our algorithm does not rely on a single global time 
step. Our algorithm avoids the requirement of an 
acute ground mesh by carefully adapting the duration 
of space-time elements to the quality of the underlying 
simplices in the space mesh. 

3. THE ADVANCING FRONT 

Our algorithm is designed as an advancing front 
procedure, which alternately constructs a patch of the 
mesh and invokes a space-time discontinuous Galerkin 
method to compute the required solution within that 
patch. To simplify the algorithm description, we 
assume that the wave speed is constant throughout 
space-time; specifically, by choosing an appropriate 
time scale, we will assume that c(p) = 1 everywhere. 
Our algorithm can be easily adapted to handle chang- 
ing wave speeds, provided the wave speed at any point 
is a non-increasing function of time. We discuss the 
necessary changes for non-constant wave speeds at the 
end of Section pi 

The input to our algorithm is a simplicial ground 
mesh M of some spatial domain X C JR d , with the 
appropriate initial conditions stored at every element. 
The advancing front M is the graph of a continuous 
time function t : X — » 1R whose restriction to any 
element of the ground mesh is linear. Any any stage 
of our algorithm, each element of the front satisfies 
the cone constraint ||Vt|| < 1. We will assume the 
initial time function is constant, but more general 
initial conditions are also permitted. 

To advance the front, our algorithm chooses a vertex 
that is a local minimum with respect to time, that 
is, a vertex p = (p,t(p)) such that t(p) < t(q) for 
every neighboring vertex q. (Initially, every vertex on 
the front is a local minimum.) To obtain the new 
front, this vertex is moved forward in time to a new 
point p" = (p,t'(p)) with t'(p) > t(p). We call the 
volume between the the old and new fronts a tent. 
The elements adjacent to p on the old front make 
up the inflow boundary of the tent; the corresponding 
elements on the new front comprise the patch's outflow 
boundary. We decompose the tent into a patch of 
simplicial elements, all containing the common edge 
pp' , and pass this patch, along with the physical 
parameters at its inflow boundary, to a DG solver. 
The solver returns the physical parameters for the 
outflow boundary, which we store for use as future 
inflow data. The solution parameters in the interior 
and inflow boundary of the tent can than be written to 




Figure 3. Pitching a series of tents over a planar triangulation. 



a file (for later analysis or visualization) and discarded. 
This advancing step is repeated until every node on the 
front passes some target time value. 

If the front has several local minima, we could apply 
any number of heuristics for choosing one; Ungor 
and Sheffer outline several possibilities pSj|. The 
correctness of our algorithm does not depend on which 
local minimum is chosen. In particular, if any vertex 
has the same time value as one of its neighbors, we 
can break the tie arbitrarily. Our implementation 
computes the mesh in phases. In each phase, we 
select a maximal independent set S of local minima 
and then lift each minimum in S, in some arbitrary 
order. This approach seems particularly amenable to 
parallelization, since the minima in S can be treated 
simultaneously by separate processors. 
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4. PITCHING JUST ONE TRIANGLE 

To complete the description of our algorithm, it 
remains only to describe how to compute the new time 
value for each vertex to be advanced, or less formally, 
how high to pitch each tent. We first consider the 
special case where the ground mesh consists of a single 
triangle. As we will show in the next section, this 
special case embodies all the difficulties of space-time 
meshing over general planar domains. 

Let p, q, r be three points in the plane. At any stage of 
our algorithm, the advancing front consists of a single 
triangle Apqr whose vertices have time coordinates 
t(p),t(q),t(r). Suppose without loss of generality that 
t(p) < t(q) < t(r) and we want to advance p forward 
in time. We must choose the new time value t'(p) 
so that the resulting triangle Ap'qf satisfies the cone 
constraint ||Vt|| < 1. 

To simplify the derivation, suppose q = (0, 0) and 
t(q) — 0. The time values t'(p) and t(r) can then 
be written as t'(p) = p ■ Vt and t(r) — r ■ X7t, where 
X7t is the gradient of the new time function. We can 
write this gradient vector as 

Vt = fj,v + vn, 

where v is the unit vector parallel to the vector r, 
and n is the unit vector orthogonal to v with sign 
chosen so that n ■ p > 0. The vector fiv is just the 
gradient of the time function restricted to segment qr, 
so jj, — t(r)/\\r\\. The cone constraint implies that 
||Vf|| = v ^i 2 + v 2 < f and therefore v < yl — . 
Thus, the cone constraint is equivalent to the following 
inequality: 

t'(p)=p-Vt 

= HP ■ v + vp ■ n 



< up ■ v + + 



p,*p ■ n 



t(r) _ V\\r\\ 2 -t(ry _ 

IMI IMI 



t(r) V¥¥~~tW 

' p-r+ v ,, „ 9 \p x r\ 



\\r 



IMI 2 



Here, p x r denotes the two-dimensional cross product 
Pi T2 — P2T\, which is just twice the signed area of 
Apqr. To simplify the notation slightly, let w p denote 
the distance from p to *qr*, and define w q and w r 
analogously: 



2|Apgr| 



2|Apgr| 



2 1 Apqr | 



||r-g|| Hp — -r|| \\q -p\ 

Then the previous inequality can be rewritten as 



t'jpX^p r I VIHI 2 -^) 2 

1 yp> - 11^112 p r + ii^ii 



(i) 



More generally, if q 7^ (0, 0) and t(q) 7^ 0, the cone 
constraint is equivalent to the following inequality. 



t'ip) < t(q) + M (p-q).(r-q) 



\\ r - all 



y/\\r-q\\*-(t(r)-t(q))\ 

+ p 



(2) 



We have similar inequalities for every other ordered 
pair of vertices, limiting how far forward in time the 
lowest vertex can be moved past the middle vertex. 
We will collectively refer to these six inequalities as 
the cone constraint. 

To ensure that our algorithm can create a mesh up 
to any desired time value, we must also maintain the 
following progress invariant: 

The lowest vertex of Apqr can always 
be lifted above the middle vertex without 
violating the cone constraint. 

This invariant holds trivially at the beginning of the 
algorithm, when t(p) = t(q) = t(r) = 0. Let us 
assume inductively that it holds at the moment we 
want to lift p. Ungor and Sheffer [T9) proved that 
if Apqr is acute, then satisfying the cone constraint 
automatically maintains this invariant, but for obtuse 
triangles, this is not enough. 

To maintain our progress invariant, it suffices to ensure 
that in the next step of the algorithm, the new lowest 
vertex q can be lifted above f without violating the 
cone constraint. In other words, if we replace t(q) with 
t(r), the new triangle's slope must be strictly less than 
f. By substituting t(r) for t(q) in the cone constraint 
(^) and making the inequality strict, we obtain the 
following: 



t'{p) < t(r) + w p 



(3) 



We have similar inequalities for every other ordered 
pair of vertices, limiting how far forward in time the 
lowest vertex can be moved past the highest vertex. 
We will collectively refer to these six inequalities as 
the weak progress constraint. 

The weak progress constraint has a simple geometric 
interpretation, which we can see by looking at the 
lifted triangle in space-time; see Figure ^. Let F be 
the cone of dependence of the lifted point p'; this cone 
intersects the plane t = t(r) in a circle 7 of radius 
t'(p) — t(r). Any plane n through p' that satisfies the 
cone constraint is disjoint from F; in particular, the 
intersection line of n with the plane t — t(r) does not 
cross 7. Now let q = (q,t(r)). If the plane p'q'f 
satisfies the cone constraint, then the line through q' 
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Figure 4. If the circle around p does not touch the line through 
q and r, then q can be lifted above f in the next step. 



and f does not cross 7. Thus, the progress invariant 
holds after we lift p only if t'(p) — t(r) < w p . 

Our algorithm lifts p to some point p' that satisfies 
both the cone constraint and the weak progress con- 
straint, where t'(p) > t(q). By the progress invariant, 
this does not violate the cone constraint. If t'(p) > 
t(r), then the weak progress constraint implies that 
the progress invariant still holds. If t'(p) < t(r), 
then the progress invariant also still holds, because 
t(r) - t'(p) < t(r) - t(q). Thus, by induction, the 
progress invariant is maintained at every step of our 
algorithm. 

Unfortunately, the weak progress constraint does not 
guarantee that we can reach any desired time value; 
in principle, the advancing front could converge to 
some finite limit. To guarantee significant progress at 
every step of the algorithm, we need a slightly stronger 
constraint. Our implementation uses the inequality 



t'{p) < t(r) + (1 - e)w p 



(4) 



where e is a fixed constant in the range < e < 1/2. 
We have a similar inequality for every other ordered 
pair of vertices, and we collectively refer to these six 
inequalities as the progress constraint. 

With this stronger constraint in place, we have the 
following result. 



Because e > and w r < ||r|| = \\r — q\\, we have 

t(rf < (l-e)V < (1-£ 2 )IM| 2 , 
which implies that 



e < 



Vlkll 2 -*M 2 



Finally, because t(r) > 0, we have 



t'(p) = ew p < 



V\\r\\ 2 -t(r) 2 



t(r) V\\r\\ 2 -t(ry 
< - ,,» p ■ r + — ^—r, 



HHP IHI 
Thus, the cone constraint is also satisfied. 



□ 



Theorem 2. Given any three points p, q,r £ 1R 2 , any 
real value T > 0, and any constant < e < 1/2, our 
algorithm generates a tetrahedral mesh of the prism 
Apqr x [0, T], where every internal facet satisfies the 
cone constraint. The number of tetrahedra is at most 
TP/2Ae, where P is the perimeter and A is the area 
of Apqr. 

Proof: Our algorithm repeatedly lifts the lowest ver- 
tex of Apqr to the largest time value satisfying the 
cone constraint (^), the progress constraint (Q), and 
a termination constraint t < T. Each time we lift a 
point, our algorithm creates a new tetrahedron. By 
Lemma jj], a new point becomes the lowest vertex, 
so the algorithm halts only when all three vertices 
reach the target plane t = T. Moreover, whenever 
t(p) < t(q) < t(r), the algorithm chooses a new 
time value t'(p) > t(q) + ew v > t(p) + ew p , except 
possibly when t'(p) — T. Thus, p is lifted at most 
T/eWp = T\\q — r\\/2Ae times before the algorithm 
terminates. □ 



Lemma 1. If the cone constraint and progress con- 
straint hold beforehand, we can lift p at least ew p 
above q without violating either constraint. 

Proof: Without loss of generality, assume that q — 
(0, 0) and t(q) = 0. We want to prove that setting 
t'(p) = ew p does not violate the cone constraint (in 
its simpler form (Q)) or the progress constraint (Q). 
Recall our assumption that t(r) > t(q) = 0. Because 
e < 1/2, we have 

t'(p) = sw p < (1 — e)w p < t(r) + (1 - e)wp, 

so the progress constraint is satisfied. The previous 
progress constraint implies that t(r) < (1 — e)w r . 



5. ARBITRARY PLANAR DOMAINS 

We now extend our meshing algorithm to more com- 
plex planar domains. The input is a triangular ground 
mesh M of some planar domain X. As we described 
in Section ^, our algorithm maintains a polyhedral 
front M with a lifted vertex p — (p,t(p)) for every 
vertex p £ M. To advance the front, our algorithm 
chooses a local minimum vertex p and lifts it to a new 
point p = (p,t'(p)). 

The new time value t'(p) is simply the largest value 
that satisfies the cone constraints and progress con- 
straints for every triangle in the ground mesh that 
contains p. The chosen time value t'(p) is the value 
that would be chosen by at least one of these triangles 
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in isolation. It follows that ft 1 is not a local minimum 
in the modified front. Moreover, by our earlier 
arguments, the progress invariant is maintained in 
every triangle adjacent to p. It follows immediately 
that our algorithm can generate meshes to any desired 
time value. 

Specifically, let u) v denote the minimum distance 
from p to *qr, over all triangles Apqr in the ground 
mesh. Lemma [l] implies the following result. 

Theorem 3. Given any triangular mesh M over any 
domain X C IR 2 , any real value T > 0, and any 
constant < e < 1/2, our algorithm generates a space- 
time mesh for the domain X x [0, T]. The number of 
patches is at most (T/e)^2 eM 1/lo p , and number of 
tetrahedra is at most (6T/e) ^2 p£M l/u P - 

Our analysis of the number of patches and elements 
is conservative, since it assumes that each step of the 
algorithm advances a vertex by the minimum amount 
guaranteed by Lemma We expect most advances to 
be larger in practice, especially in areas of the ground 
mesh without large angles. Our experiments were 
consistent with this intuition; see Section |^. 

Most of the parameters of the cone constraint, and 
all of the parameters of the progress constraint, can 
be computed in advance from the ground mesh alone. 
Thus, the time to compute each new time value t'(p) 
is a small constant times the degree of p in the ground 
mesh, and the overall time required to build the mesh 
is a small constant times the number of mesh elements. 

Non-constant Wave Speeds 

Although we have described our algorithm under the 
assumption that the wave function c(p) is constant, 
this assumption is not necessary. If elements of the 
ground mesh have different (but still constant) wave 
speeds, our algorithm requires only trivial modifi- 
cations. The situation fits well with discontinuous 
Galerkin methods, which compute solutions with dis- 
continuities at element boundaries. If the wave speed 
varies within a single element, even discontinuously, 
the only necessary modification is to compute and use 
the maximum wave speed over each entire element. 
Similar modifications suffice if the wave speed at any 
point in space can decrease over time. 

If the mesh has only acute angles, the progress 
constraint is redundant and arguments of Ungor and 
Sheffer [Jl9| imply that our algorithm works even if the 
wave speed can increase over time, as long as the wave 
speed is Lipschitz continuous. Unfortunately, their 
analysis breaks down for obtuse meshes because of the 
progress constraint, and indeed our algorithm can get 



stuck. We expect that a refinement of our progress 
constraint would allow for increasing wave speeds, but 
further study is required. 

6. HIGHER DIMENSIONS 

Our meshing algorithm extends in an inductive man- 
ner to simplicial meshes in higher dimensions. As in 
the two-dimensional case, it suffices to consider the 
case where the ground mesh consists of a single simplex 
A in IR d . At each step of our algorithm, we increase 
the time value of the lowest of the simplex's d + 1 
vertices as much as possible so that the cone constraint 
||Vi| < 1 is satisfied and we can continue inductively 
as far into the future as we like. 

Let p, 5, n, J"2, . . . , r<j-i denote the vertices of A in 
increasing time order, breaking ties arbitrarily. Our 
goal is to lift p above q without violating the cone 
constraint. Let F be the facet of A that excludes p, 
and let H be the hyperplane spanning F. Let p H be 
the projection of p onto H , and let p F be the closest 
point in F to p. Observe that Zpp H p F is a right angle. 
See Figure ^. Let V H t denote the gradient vector of 
the time function restricted to H . Finally, define 

\\p-Ph\\ 
Wv-pAY 

If p H lies inside F, then p H = p F and a F — 1; otherwise, 
p H / p F and a F = sin Zp H p F p. 



P 




Figure 5. Defining the points p H , p F , and p z 



The higher-dimensional analogue of the weak progress 
constraint is described by the following lemma. 

Lemma 4. If \\V H t\\ < a F , then we can lift p above q 
without violating the cone constraint ||V£|| < 1. 

Proof: Suppose ||V^i|| < o F . Without loss of gener- 
ality, assume that q = (0, 0, . . . , 0) and t(q) = 0. Let 
n be the unit normal vector of H with p • n > 0, 

V -P H 
n = li ir 

Since the time function t is linear, changing only t(p) is 
equivalent to leaving t fixed on the hyperplane H and 
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changing the directional derivative dt/dn. To prove 
the lemma, we show that setting 

|i = cosZ OT =M^M (5) 
dn " * lb -fell 

gives us a new time function that satisfies the cone 
constraint with t(p) > 0. 

Let Z be the set of points in H where t = 0. Since 
t(q) = 0, Z is the (d — 2)-flat orthogonal to V H t that 
passes through q. Moreover, because t > everywhere 
in F, Z is a supporting (d — 2)-flat of F. Let p z be the 
closest point in Z to p (or to p H ); this might be the 
same point as p F , p H , or q. Observe that Zpp H p z is a 
right angle. See Figure 

We can express the time gradient Vt as follows: 

^ ^ 9t 
Vt = V„£ + — n. 

an 



Equation (JE]) implies that 



Vi = V w i + 



llftf-Pf 



Since these two components of are orthogonal, we 
can express its length as follows. 

l|Vt|| 2 = ||V„i|| 2 + ^^£ 



IIP ~Pf 



IIp-PhII : 



+ 



= l 



IIp-PfII 2 Hp -fell 

So the new time function satisfies the cone constraint. 
We can express the time value t(p) as follows: 

t(p) = -t(p H ) + \\p - Ph\\ 

= t(p H ) + 



dt_ 

dn 

p H \\ \\p h - 
lb-fell 



If t(p H ) > 0, then clearly t(p) > 0. Suppose t(p H ) < 0. 
The vector p H — p z is orthogonal to Z and therefore 
anti-parallel to V H t. Thus, 

t(p H ) =V H t- (p H -P z ) 

= -||V H *|| life -fe|| 

^ lb -fell life -fell 



Hp-PpII 
■fell life -fell 



lb -fell 

The last inequality follows from the fact that p H and F 
lie on opposite sides of Z, because t(p H ) < 0. It now 
immediately follows that t(p) > 0. □ 

As in the two-dimensional case, in order to guarantee 
that the algorithm does not converge prematurely, 
we must strengthen this constraint. There are many 
effective ways to do this; the following lemma describes 
one such method. 



Lemma 5. For any < £ < 1, if ||V H £|| < (1 — s)cr F , 
then we can lift p at least e\\p — p H \\ above q without 
violating the cone constraint ||Vt|| < 1. 

Proof: We modify the previous proof as follows. We 
show that setting 

dn +l J Hp-p f II 

gives us a new time function satisfying the conditions 
of the lemma. First we verify that the cone constraint 
is satisfied. 

m f < ( (1 _ £ )t^H) 2 + (s + (1 - 



i + 2e (i- £ )fM^M-i 
b-fe 



< i 



(In fact, if p H 7^ p F , then ||Vt|| < 1, which means we 
could lift p even more.) 

Next we verify that t(p) > e\\p — p H \\. 

at 



t(p)=t(p H ) + \\p-p H \\gr 

= t(p H ) + e\\ P -p H \\ + (l-e) 
>t(p H )+e\\p-p H \\. 



IIp-PhII IIph-PfII 
lb -fell 



If t(p H ) > 0, we are done. Otherwise, as in the previous 
lemma, we have 

«(&)>- II V H t|| ||&-fe|| 

lb -fell life -fell 



>-(i-^) J 



Hp - Pf 



which immediately implies that t(p) > e\\p — p H \\, as 
claimed. □ 



An important insight is that we can view the simplex 
A simultaneously as a single d-dimensional simplex 
and as (d — l)-dfmensional boundary mesh. Lemma [B| 
prescribes a tighter cone constraint for every element 
in this boundary mesh. 

Our algorithm proceeds as follows. At each step, we 
lift the lowest vertex of A by recursively applying the 
(d — l)-dimensional algorithm; then, if necessary, we 
lower the newly-lifted vertex to satisfy the global cone 
constraint ||V£|| < 1. The base case of the dimensional 
recursion is the two-dimensional algorithm in the 
previous section. 

This recursion imposes an upper bound on the length 
of the time gradient within every face of A of dimen- 
sion at least 1. In fact, a naive recursive implementa- 
tion would calculate (d — k)\ different constraints for 
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each fc-dimensional face. A more careful implementa- 
tion would determine the strictest constraint for each 
face in an initialization phase, so that each step of the 
algorithm only needs to consider each face incident to 
the lifted vertex once. 

For a d-dimensional ground mesh with more than one 
simplex, we apply precisely the same strategy as in the 
two-dimensional case. At each step of the algorithm, 
we choose an arbitrary local minimum vertex p, and 
lift it to the highest time point p allowed by all 
the simplices (of all dimensions) containing p. By 
our earlier arguments, p' is not a local minimum of 
the modified front, which implies that our algorithm 
terminates only when all the vertices reach the target 
time value. 

7. OUTPUT EXAMPLES 

We have implemented our planar space-time meshing 
algorithm and tested it on several different ground 
meshes. Our implementation consists of approxi- 
mately 5000 lines of C++ code, about 800 of which 
represent the actual space-time meshing algorithm; 
the remaining code is a pre-existing library for ma- 
nipulating and visualizing triangular and tetrahedral 
meshes. 

Figures ^ and ^|-^| show space-time meshes computed 
by our implementation. In each case, we stopped 
advancing each vertex of the front after it passed a 
target time value. In every example, the input triangle 
mesh contains at least one (sometimes extremely) 
obtuse triangle, which caused Ungor and Sheffer's 
original Tent Pitcher algorithm to fail |19| . 

Our program produces several thousand elements per 
second, running on a 1.7 GHz Pentium IV with 
1 gigabyte of memory. For example, the mesh in 
Figure hi, which contains 114,515 tetrahedral elements, 
was built from a ground mesh of 2,356 triangles in 
about 14 seconds. Figure ^j] shows an input mesh with 
1,044 triangles and the resulting 55,020-element space- 
time mesh, which was computed in about 4 seconds. 
(These running times include reading and parsing the 
input mesh file and writing the output mesh to disk.) 

Figure ^ illustrates effect of grading in the input mesh 
on the size on space-time elements. The largest and 
smallest elements in the ground mesh differ in size by 
a factor of 128; the resulting space-time elements differ 
in duration by a factor of 450. (The difference between 
these two factors might be explained by the obtuse 
triangles near the smallest element of the ground 
mesh.) Less severe grading due to varying ground 
element size can also be seen in Figure ^. 

Figure || shows the output of our algorithm when the 
input mesh is pathological. The input meshes are 




(b) 

Figure 6. (a) A typical planar mesh of 1044 triangles, (b) The 
resulting space-time mesh of 55,020 tetrahedra, computed by 
our implementation in about 4 seconds. 



the Delaunay triangulation and a greedy sweep-line 
triangulation of the same point set. As expected, 
variations in quality in the ground mesh also leads to 
temporal grading in our output meshes. For example, 
the bottom right vertex of the space mesh in Figure 
^|(b) advances much more quickly than the top right 
vertex, because it is significantly further from the lines 
through any of its neighboring edges. 

We tried several different values of the parameter e 
in the progress constraint (Q). All of the example 
output meshes were computed using the value e « 0.1. 
Somewhat to our surprise, the number of elements 
in the output mesh varied by only a few percent as 
we varied e from 1/100 to 1/3, and smaller values 
of e usually resulted in meshes with slightly fewer 
elements, since the modified progress constraint is less 
severe. Also, for high-quality ground meshes, where 
most of the triangles are acute, the progress constraint 
affected only a few isolated portions of the space- 
time mesh. On the other hand, smaller values of e 
generally led to wider variability in the duration of 
neighboring tetrahedra. As e increases, the progress 
guaranteed by Lemma [j] more closely matches the 
maximum progress allowed by the progress constraint; 
this tends to distribute the progress of each triangle 
more evenly among its vertices. 
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(c) 

Figure 7. (a) A severely graded planar mesh, (b) The resulting 
space-time mesh, (c) A close-up of the resulting grading. 



8. FURTHER RESEARCH 

We have presented the first algorithm to generate 
graded space-time meshes for arbitrary spatial do- 
mains, suitable for efficient use by space-time dis- 
continuous Galerkin methods. This is only the first 
step toward building a general space-time DG meshing 
library. 

As we mentioned in Section |j| our algorithm currently 
requires the wave speed at any point in space to remain 
constant or monotonically decrease over time. In the 
short term, we plan to adapt our algorithm to handle 
wave speeds that increase over time. It should be 
noted that for many problems, the wave speed is not 
known in advance but must be computed on the fly as 
part of the numerical solution. 

DG methods do not require conforming meshes, where 
any pair of adjacent elements meet in a common face. 
As a result, fixed time-step methods allow the space 
mesh to be refined or coarsened in response to error 
estimates, simply by remeshing at any time slice. 
Can our advancing front method be modified to allow 
for refinement, coarsening, or other local remeshing 
operations (like Delaunay flips)? These operations 
might be useful not only to avoid numerical error, but 
also to make the meshing process itself more efficient. 

For many problems, even the boundary of the domain 
changes over time according to the underlying system 
of PDEs. Can our method be adapted to handle 
moving boundaries? Intuitively, we would like a mesh 
that conforms to the boundary as it moves. This 
would require us to move the nodes of the ground mesh 
continuously over time; remeshing operations would be 
required to guarantee that the meshing algorithm does 
not get stuck. Similar issues arise in tracking shocks, 
which are surfaces in space-time where the solution 
changes discontinuously. 

Our method currently assumes that all the charac- 
teristic cone have vertical (or at least parallel) axes. 
For problems involving fluid flow, the direction of the 
cone axis (i.e., the velocity of the material) varies 
over space-time as part of the solution. We could 
adapt our method to this setting by overestimating 
the true tilted influence cones by larger parallel cones, 
but intuitively it seems more efficient to move nodes 
on the front. As in the case of moving boundaries, this 
would require remeshing the front. In fact, the front 
would no longer necessarily be a monotone polyhedral 
surface; extra work may be required to ensure that the 
resulting mesh is acyclic. 

Finally, to minimize numerical error it is important to 
generate space-time meshes of high quality. Although 
there are several possible measures for the quality of 
a space-time element, further mathematical analysis 
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(c) (d) 

Figure 8. (a) A Delaunay triangulation with a few bad triangles, (b) A sweep-line triangulation (of the same point set) with many 
horrible triangles. (c,d) The resulting space-time meshes, showing the resulting temporal grading. 



of space-time DG methods is required to determine 
the most useful quality measures. This is in stark 
contrast to the traditional setting, where appropriate 
measures of quality and algorithms to compute high- 
quality meshes are well known jlj ^, |l3|, [L5| . 
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